Violation of the fluctuation-dissipation theorem in a protein system 
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We report the results of molecular dynamics simulations of the protein myosin carried out with 
an elastic network model. Quenching the system, we observe glassy behavior of a density correla- 
tion function and a density response function that are often investigated in structure glasses and 
spin glasses. In the equilibrium, the fluctuation-response relation, a representative relation of the 
fluctuation-dissipation theorem, holds that the ratio of the density correlation function to the density 
response function is equal to the temperature of the environment. We show that in the quenched 
system that we study, this relation can be violated. In the case that this relation does not hold, this 
ratio can be regarded as an effective temperature. We find that this effective temperature of myosin 
is higher than the temperature of the environment. We discuss the relation between this effective 
temperature and energy transduction that occurs after ATP hydrolysis in the myosin molecule. 

PACS numbers: 87.14.Ee, 05.70.Ln, 05.40.-a 



I. INTRODUCTION 

Proteins possess complex structures and, consequently, 
complex motion. Such complexity might be necessary to 
carry out the functions exhibited by living organisms. 
Statistical properties of the types of complex structure 
and motion characterizing protein molecules have been 
studied using various approaches, in particular, that em- 
ploying energy landscapes [lj. Recently, with a compu- 
tation of the density of states for a Go-like model of a 
protein, statistical properties were investigated using the 
idea of inherent structures (2j. An inherent structure is 
a subset of the configuration space that represents the 
local minima of an energy landscape, originally proposed 
to study the dynamics of liquids [3j . Although the idea 
of inherent structures is very important to understand 
protein dynamics, it is difficult to directly experimen- 
tally investigate the inherent structures of a protein. For 
this reason, it would be useful if we could characterize 




FIG. 1: Schematic depiction of the myosin molecule (1KK7 
[I?]])- The horizontal axis denotes the x direction, and the 
vertical axis denotes the y direction. According to our con- 
vention, the a-carbon atoms, C a with with i = 1, • ■ • , 770 
belong to the head substructure (white), and those with 
i — 771, • ■ ■ , 1101 belong to the tail substructure (gray). See 
the next section for the way how to number C a atoms. 



the energy landscape of a protein in forms of experimen- 
tally measurable quantities, such as a density correlation 
function or a density response function. 

In structure glasses, which, like protein molecules, also 
possess energy landscapes with many local minima, sta- 
tistical properties are often studied by computing the 
density correlation function and the density response 
function 0, i, i, 0, 0, H- In equilibrium, these quan- 
tities satisfy the fluctuation-response relation, a repre- 
sentative relation in the fluctuation-dissipation theorem 
10]. This relation means the equality of the ratio of the 
density correlation function to the density response func- 
tion and the temperature of the environment [10]. For 
glassy systems, which exhibit slow relaxation and are in- 
herently non-equilibrium, it has been reported that the 
fluctuation-response relation is violated 0, H, H, @, H, [1] ; 
that is, the ratio of the relevant density correlation func- 
tion and density response function is not equal to the 
temperature of the environment. The slow relaxation 
displayed by a glassy system, which results from the na- 
ture of its energy landscape, with many local minima, 
can be characterized by a quantity representing the de- 
gree of violation of the fluctuation-response relation. In 
some cases that the fluctuation-response relation is vio- 
lated, the ratio of the density correlation function to the 
density response function has interpreted as an "effec- 
tive temperature," and there are studies addressing the 
question of whether this effective temperature can play 
the role of the temperature in non-equilibrium systems 

In this paper, we report the results of simulations of 
molecular dynamics employing an elastic network model 
of sub-fragment 1 of a myosin molecule, which is com- 
posed of a head substructure, with ATP-binding and 
actin-interacting sites, and a tail substructure, with a 
long alpha- helix bound by two light chains (see Fig. [1]). 
With the hypothesis that appropriate density correlation 
function and density response function, which have not 
yet been fully exploited in studies of protein dynamics, 
can be used to characterize the glassy behavior of pro- 



2 



tein molecules, we compute such quantities in the case of 
a quenched system. While elastic network models have 
been used to study the elastic properties of equilibrium 
fluctuations [111 Q2, such models have not been used 
to study non-equilibrium behavior. In this paper, we 
show that glassy, non-equilibrium behavior can also be 
described with this class of models. 

Because the effective temperature that we employ rep- 
resents the degree to which the behavior of the system 
is "glassy," we quantitatively investigate the complexity 
of the myosin's head through use of this effective tem- 
perature. More specifically, we compute the values of 
the effective temperature for both the head and the tail, 
and we compare these values with the temperature of 
the environment. We also seek to elucidate the origin 
of the observed glassy behavior through analysis that in- 
vestigates the inherent structures We find that the 
inherent structure for the model we study has a resem- 
blance to the structure of a structural isomer of myosin 
determined by X-ray crystallography. 



Protein Data Bank [l6j] ■ Here, we chose the values of rf 
so as to obtain the structure of the myosin molecule I KK7 
[13] • Also, these values are such that the center of mass 
is at (0, 0, 0) and the unit of length is A. 

The function Vt va , v ({ri}) that appears in |T]) is a trap- 
ping potential that plays the role of an optical potential, 
acting to fix the molecule. This potential is given by 
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The time evolution of the system is described by the 
Langevin equation (i = 1, • • • , TV) 



dvi dH . . 
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dt 
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where £j = (£«, £ Z i) is Gaussian white noise that 
satisfies 



II. ELASTIC NETWORK MODEL 

In our study, a protein molecule is regarded as consist- 
ing of a-carbon atoms, C Q , employin g a coarse-grained 
representation of amino acid residues [13L IbH . [lj| . A C a 
atom is a constituent of the carbon skeleton of a protein 
molecule. 

The position and velocity of the i-th C a in a myosin 
molecule are denoted by = (xi,yi,Zi) and Vi = 
(v X i,v y i,v Z i), where i = 1, ••• ,N, and N is the total 
number of C Q atoms in the myosin molecule. Here, C Q 
atoms are numbered in the order of the heavy chain, the 
essential light chain and the regulatory light chain. 

The Hamiltonian of our model is given by 

JV 

H({r l }AM) = Y, ! v\ v ^ + V ({r l }) + V tTap ({r l }), (1) 



where m is the mass of the i-th C a , and V({ri}) is the 
interaction potential 
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We stipulate that k e = (£ = 1,2,3) if |r° - r°| > r c , 
where r c is a cut-off length, and we set = 0.5fci and 
&3 = O.lfci. With this form of V({ri}), the configuration 
of the native structure, {r"}, is most stable. The values 
of the parameters r° used here were taken from the RCSB 



(O) = 2 7^(i - t')6 a ,fi5 hJ 



(6) 



Here, T is the temperature of the environment (with the 
Boltzmann constant set to 1), and 7 is the friction con- 
stant of the solvent. Here, ( ) denotes the average over all 
noise histories. We set the parameters used in our numer- 
ical simulations as m = 1, 7 = 0.01, ki = 1, N = 1101 
and r c = 10. 



III. EQUILIBRIUM BEHAVIOR: 
FLUCTUATION-RESPONSE RELATION AT 
HIGH TEMPERATURE 

In equilibrium, the ratio of the density correlation 
function to the density response function is equal to 
the temperature of the environment. This relation 
among the density correlation function, the density re- 
sponse function and the temperature of the environ- 
ment is called the fluctuation-response relation (For a 
review of the fluctuation-dissipation theorem, including 
the fluctuation-response relation, see Ref. jlOj). 

In order to examine whether the fluctuation-response 
relation holds for our system, we first introduce a density 
response function. In this section, we consider the relax- 
ation process that results in the case that the system is 
initially in equilibrium, and then at t — the perturbing 
potential 



N 



N 



E y p(^) = E Acos ^) 



(7) 



(=1 



is added to the Hamiltonian (fT]). Note that we consider 
a perturbation along the y direction, which is approxi- 
mately parallel to the long axis of the myosin (see Fig. 
[I}. We have also studied the cases with perturbations 
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Then, the following is one representation of 
fluctuation-response relation: 



the 



R(t) = 



1 dC(t) 
T dt 



(11) 



where t > [13]. Here, the density response function, 
R(t), defined as the time derivative of x(i): 



R(t) 



d X (t) 
dt ' 



(12) 



In numerical experiments, we computed x(t) m order to 
investigate the behavior of R(t). 

In the upper graph of Fig. [2j xW an d C(t) are plot- 
ted as functions of time with T = 0.8, k = 27r/10 and 
A = 0.05. In the lower graph of Fig. [2J xif) is plotted as 
a function of C(t). The fact that the slope of the lower 
graph is — 1/T indicates that the fluctuation-response re- 
lation (jlip holds in equilibrium (the high temperature 
regime). Note that we set k — 27r/10 here because the 
lengths of the alpha helices in myosin are of the order of 
10 A. Glassy behavior of the density correlation and the 
density response is observed in the case k = 2tt/£ with 
I > 10, as described below. 



FIG. 2: (Top) C(t) (circles) and x(*) (triangles) plotted as 
functions of time in the equilibrium case with T — 0.8, A = 
0.05 and k = 2-7r/10. (Bottom) x(t) as a function of C(t) 
in the case T = 0.8. The slope of the line here is —1/T. 
In both graphs, the data points are obtained by averaging 
2000 independent trajectories, and the statistical error bars 
are smaller than the symbols. 



along different directions, and the results obtained in all 
cases are qualitatively the same as those obtained along 
the y direction. In this relaxation process, the suscepti- 
bility, x(*)> is defined as 



x(t) 



(p(M)-p(M)r p 



(8) 



where ( } p denotes the statistical average under this re- 
laxation process. Writing the density in the y direction 

as p(y,t) = YliLi^iv ~ Vi(t))/N, we denote its Fourier 
transform by p(k, t): 

P (k,t) EE 



dyp(y,t) cos(ky) 



1 N 

— J2cos(ky,(t)). 



(9) 



Next, we define the density correlation function, which 
is computed in equilibrium without adding V v , as 



C(t) = (p(k, t)p(k, 0)} N. 



(10) 



IV. NON-EQUILIBRIUM BEHAVIOR I: 
GLASSY PHENOMENON AT LOW 
TEMPERATURE 

In this section, we consider simulations in which we 
quenched the system from T = 0.5 to T = 0.05 in order 
to investigate its glassy behavior. In this model, room 
temperature roughly corresponds to T = 0.1, which is 
estimated by comparing the mean square fluctuations of 
the C a atoms in the model with those obtained in an all- 
atom model (data not shown). Therefore, the quenched 
temperature T = 0.05 roughly corresponds to 150 K, 
which is below the glass transition temperature of a real 
protein [181 ]. Note that myosin remains in a folded form 
even under conditions corresponding to T = 0.5. 

In our simulations, we quench the system at t = —t w , 
and then we begin investigating the behavior of the sys- 
tem at t = 0. The quantity t w is called the "waiting 
time." In general, relaxational glassy systems do not 
possess invariance with respect to time translation, and 
both correlations and responses decay more slowly as the 
age of the system increases. Thus, in such systems, cor- 
relation functions and response functions depend on t w . 

In order to elucidate the glassy behavior of our system, 
we first introduce the auto-correlation function 

F(t,t w ) ee ^1 J2[cos(y % (t) - W (0)) - cos(y° - w (0))]\ , 

(13)*™ 

where ( ) t represents the statistical average under the 
relaxation we consider. (Auto-correlation functions are 
often employed in the study of glassy systems 0, [H, H, 0, 
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FIG. 3: F(t,t w ) as a function of time for the case i w = 100. 
The data points are obtained by averaging 2000 independent 
trajectories, and the statistical error bars are smaller than the 
symbols. 



& Q ■) The quantity F(t, t w ) is a measure of the similarity 
of the configuration of the system at time t to that at time 
0; note that the quantity cos (yi(t) — j/i(0)) is equal to 1 
at t — and decreases as t increases. The second term on 
the righthand side of (|13|) ensures that F(t,t w ) = for 
large t, because the configuration of the system, {^(i)}, 
fluctuates about that of the native structure, {r°}, at 
large t. Computing F(t,t w ), we canelucidate the slow 
relaxation of the configuration of the system. In Fig. [3J 
we plot F(t,t w ) as a function of time for the case f w = 
100. The form of the decay, which is not exponential, 
reveals the slow relaxation of the configuration. 



V. NON-EQUILIBRIUM BEHAVIOR II: 
VIOLATION OF THE 
FLUCTUATION-RESPONSE RELATION AT 
LOW TEMPERATURE 

In an aging system, the susceptibility, density response 
function and density correlation function, which depend 
on t w , are defined by 



x(Mw) = - 



(p(M)-p(*,o))£ 



R{t, t v 
C(t,t y 



dx(Mw) 
dt ' 
(p(k,t)p(k,0)). n. 



(14) 

(15) 
(16) 



In the lower graph of Fig. [H we plot x(Mw) as a 
function of C(t, t w ) for the cases i w = 100 and i w = 200. 
It is seen that there are two slopes, characterizing two 
time regimes, unlike in the equilibrium case (see the lower 
graph in Fig. [2]). The results plotted in Fig. [2] suggest 
that in the early regime, the form of x(t, t w ) as a function 
of C(t,t w ) approaches a line of the slop — 1/T, where 
T = 0.05, as t w increases. (See Fig. 16 of Ref. [§] for 
such observation of similar behavior.) In the late time 
regime, the form of this function is again a line, but in 
this case, the slope differs significantly from —1/T. (Also 
note that in this case, similar linear behavior is observed 
for both < w = 100 and t w — 200.) Such behavior of 
x{t,t w ) and C(t,t w ) is often observed in aging systems, 
including structure glasses 0, [f| and spin glasses [J, [f| . 
In an aging system, the slope of \(t, t w ) as a function 
of C(t,t w ) in the late regime is — 1/T e ff where T e g- is 
called the effective temperature, and is defined through 
the relation 



R(t, i w ) 



1 dC(*,t w ) 



T, 



cir 



dt 



(17) 



We interpret this behavior as indicating that in the early 
regime, the system relaxes toward the local equilibrium 
represented by a meta-stable state, while in the late 
regime, the system wanders among many meta-stable 
states, as it evolves toward the minimum of the free en- 
ergy of the system. These two kinds of relaxational pro- 
cesses exhibited by an aging system are characterized by 
the two slopes —1/T and — l/T e g. 

The fact that T c g ^ T represents a violation of the 
fluctuation-response relation in this aging system. The 
lack of invariance with respect to time translation and 
the violation of the fluctuation-response relation are a 
basic properties of relaxational glassy systems. The effec- 
tive temperature T c g denotes the difference of the system 
from equilibrium. 

Before ending this section, we note an important point 
regarding T e g. The inequality TsB> T is observed in 
many glassy systems 0, d, H, @, H II]- However, because 
there are counterexamples, the universality of this in- 
equality has not yet been established. Here, we found 
T Q g > T in the elastic network model. Possible implica- 
tions of this inequality in biological contexts are discussed 
below. 



VI. T cB FOR THE HEAD AND TAIL 



Note that ( ) t p denotes the statistical average under the 
condition that we quench the system at t = — t w , and 
switch on the perturbation potential given by ([7]) at t — 
0. In the upper graph of Fig. 5, C(t,t w ) and x(i,i w ) 
are plotted as functions of time for the case t w = 100. 
The relaxations of C(t,t w ) and x(*>*w) seem to be slow 
in comparison with those of C(t) and x(t) found in the 
equilibrium case with T — 0.8 (see the upper graph in 
Fig. 2). 



One of the biggest differences between a structure glass 
and a protein is that a protein has characteristic sub- 
structures that are believed to cause a smoothing of the 
energy landscape, which makes the collective motion nec- 
essary for biological functions possible. As mentioned 
above, a myosin molecule is composed of a head, with 
ATP-binding and actin-interacting sites, and a tail, with 
a long alpha-helix bound by two light chains. It is inter- 
esting to compare these two parts with respect to T e ff. 
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FIG. 4: (Top) C(Mw) (circles) and x(Mw) (triangles) plotted 
as functions of time in the case of an aging system, with t w = 
100, and A = 0.05 and k = 2tt/10. (Bottom) x(M») as 
a function of C(t,t w ) for the cases t w = 100 (circles) and 
t w = 200 (triangles). The slope of the dotted line is —1/T, 
where T = 0.05. The slope of the solid line is — 1/T e s, where 
T c s = 1.3. In both graphs, the data points are obtained by 
averaging 2000 independent trajectories, and the statistical 
error bars are smaller than the symbols. 

For this purpose, we investigated the density response 
function and density correlation function for the head 
(i = 1, ■ ■ ■ , 770) and the tail (i = 771, • • • ,N), and com- 
puted T e ff for each. 

We define the Fourier transform of the density for the 
head and tail as 

770 

Ph{k,t) ee — ^cos(fcy l ), (18) 

l — l 

1 - 

Pt(k,t) ee N _ 77Q C0S ( fc ^)- (19) 

i=771 

Then, we use p\(k,t) and p t (k,t) instead of p(k,t) in 
dHJ) and (JUS). Also, letting N h = 770 and N t = N — 
770 be the total numbers of Cq, atoms in the head and 
tail, respectively, we use these in place of N in ([T3|) and 
(fTB]) . We note here that when computing the response 
functions, we add a perturbing potential of the form given 
in (J7J only to the substructure under investigation. 
In Fig. we plot x(Mw) as a function of C(t,t w ) 
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FIG. 5: %(£, t w ) as a function of C(t, t w ) in the case t„ = 100 
and k = 27r/10 for the head (circles) and the tail (triangles) 
of the myosin molecule. The slope of the solid line is equal to 
— 1/1.3, while the slope of the dotted line is equal to —1/0.6. 
The data points are obtained by averaging 2000 independent 
trajectories, and the statistical error bars are smaller than the 
symbols. 

for both the head and the tail, respectively. We find that 
the effective temperature, T e g , for the head is higher than 
that for the tail. This implies that the structure of the 
head is more glassy than that of the tail. 

VII. DISCUSSION 

In conclusion, we investigated the glassy behavior of an 
elastic network model of the myosin molecule by studying 
the density correlation function and the density response 
function in the case that we quench the system from T = 
0.5 to T — 0.05. The glassy behavior is displayed in Figs. 
2] and OB where the susceptibility, x(t, £ w ), is plotted as a 
function of the correlation function, C(t,t w ). We found 
that T cS defined in (p2]) was not equal to T, and thus 
that the fluctuation-response relation is violated. We also 
compare the degrees of the violation of the fluctuation- 
response relation for both the head and tail substructures 
of myosin, individually, and we found that T c g is higher 
for the head than for the tail. In the following, we discuss 
two points related to these main results. 



A. Effective temperature of the myosin molecule 

Although further studies are required to ascertain 
physically clear interpretations of the effective temper- 
ature, the fact that T e fj is higher than T is intriguing, 
because it may be related to energy transduction taking 
place after ATP hydrolysis in a myosin molecule. 

An acto-myosin system, of which a myosin molecule 
is a component, is a representative motor-protein sys- 
tem, and it has been studied in single-molecule experi- 
ments [l9l . [20l [2ll | . Upon binding of ATP to a myosin 
molecule, the myosin molecule detaches from an actin 



filament. ATP hydrolysis takes place in the detached 
myosin molecule. The products of the ATP hydroly- 
sis (ADP and an inorganic phosphate) are then released 
from the myosin molecule, which are thought to be cou- 
pled with the force exertion and with the re-attachmcnt 
of myosin to the actin filament. Recently, a single- 
molecule experiment on an acto-myosin system demon- 
strated that the force exertion of the myosin molecule 
sometimes occurs after the release of the bound ADP [lj| . 
This fact leads us to believe that the energy provided 
by the ATP hydrolysis might be stored in the myosin 
molecule for a short time before the force is exerted. If 
this is indeed the case, it is important to determine the 
form in which this energy is stored. The effective tem- 
perature, T e ff, might help us solve this problem. Because 
ATP hydrolysis takes place only in the head, and because 
T e g is higher in the head than in the tail, it is reasonable 
to conjecture that this energy storage has a close con- 
nection with the inequality T c g > T. To elucidate this 
connection is a future problem. 



B. Inherent structures 

Next, we discuss the origin of the glassy behavior 
observed in the elastic network model we study. In 
Fig. [(JJ we plot E\s, which is the total potential energy, 
V + Vtrap, when the system is at a local minimum in the 
energy landscape, as a function of time. The quantity 
E\s is computed by using the steepest-descent energy- 
minimization method, as has been employed in Ref. [3j, 
in the case t w = 100. From the rather discretized forms 
of the trajectories, it is seen that the system moves from 
one meta-stable state to another. We stress that the 
clastic network model described by @ is not a harmonic 
system, and hence can possess local energy minima, al- 
though the model may appear to be a harmonic system 
(Actually, cxistense of local energy minima in an elastic 
network model is indicated in Ref. [22j|). 

It would be interesting to examine the difference be- 
tween the native structure and the structures correspond- 
ing to the local energy minima (i.e., inherent struc- 
tures). The upper graph of Fig. [JJ displays the difference 
\yj s — yi \ between the positions of the C a atoms in the in- 
herent structure and in the structure, 1KK7, where {y] s } 
is the y value of the position of the i-th C Q atom in the 
inherent structure. From this graph, it is seen that the 
differences for some of the C Q atoms are as much as 10 
Aand that those C Q which exhibit such large displace- 
ments are not localized but, rather, clustered. 

Several structural isomers of myosin have been experi- 
mentally found by using X-ray crystallography [23] . The 
structure we studied here, 1KK7, is one with no nu- 
cleotide bound. Therefore, it would be interesting to 
compare the inherent structures with those of myosin 
isomers. As an example, we compare one of these in- 
herent structures with the structure of the myosin iso- 
mer with an ATP-analog bound, 1KK8 [l7j]. The lower 
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FIG. 6: (Top) Two examples of Eis, computed using 
the steepest-descent energy-minimization method, plotted as 
functions of time in the case t„ = 100. (Bottom) Three ex- 
amples of the steepest-descent energy-minimization trajecto- 
ries are shown. Note that the three energy-minimization tra- 
jectories, with different initial configuratio ns (i.e., different 
"instantaneous structures" on the MD trajectory), tend to 
converge at different energy values. 



graph in Fig. [7J shows that there is a resemblance be- 
tween the inherent structure and the structure of the 
myosin isomer, with relatively large displacements seen 
near the N-terminal domain, the lower-50k domain, and 
the converter domain. This suggests that the elastic net- 
work model contain information concerning the struc- 
tures of other isomers. Using a normal mode analysis 
for an elastic network model of a myosin molecule, Zhen 
and Doniach have shown that a structural isomer is lo- 
cated along the directions of some of the slowest modes 
of the structure [ll| . Our results appear to be consistent 
with their results. Furthermore, our results indicate the 
meta-stability of the isomer. We believe that our results 
may lead to an extension of the applicability of the elas- 
tic network model of proteins, noting that it has been 
shown that this model can even be applied to the study 
of protein folding [24| and the investigation of nonlin- 
ear relaxation dynamics [22| . To elucidate the range of 
applicability of the elastic network model, we need to sys- 
tematically study its meta-stable states for many kinds 
of proteins. 
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FIG. 7: (Top) Three examples (black, pink, green) of \y] — yf\ plotted as functions of i (1 < i < 770 , the head substructure 
in our definition), where {yj s } is a configuration computed using the steepest-descent energy- minimization method, and {y®} 
is the configuration of the structure 1KK7. Note that the a-carbon atoms, C Q , with i — 350, ■ ■ ■ , 400 are trapped by Kra P (see 
p|). (Bottom, Left) An example of \yj s — y®\ (black) and \yl KKS — yfl (red) plotted as functions of i. Here, {yl KKS } is the 
configuration of the structure 1KK8. (Bottom, Right) \yl KKa — yf \ plotted as a function of \y] s — y1\. 



